# This file replicates analysis for: 
# "How does community policing affect police attitudes? An experimental test and a theory of bureaucrat-citizen contact"
# Dotan Haim, Matthew Nanes, and Nico Ravanilla
# American Political Science Review

################################################
############### 1. Preparatory #################
################################################

	# =======================================================
	# Set up Workspace
	# =======================================================	
		# Clear Workspace
		rm(list=ls())
		
		# Install Packages
			# options(repos = c(CRAN = "https://cran.rstudio.org"))		
			# install.packages("tidyr")
			# install.packages("data.table")
			# install.packages("estimatr")
			# install.packages("ggplot2")
			# install.packages("gridExtra")
			# install.packages("texreg")
			# install.packages("kableExtra")
			# install.packages("vtable")
			# install.packages("margins")
			# install.packages("purrr")
			# install.packages("interflex")

		# Load Packages
		suppressMessages(library("tidyr"))
		suppressMessages(library("data.table"))
		suppressMessages(library("estimatr"))
		suppressMessages(library("ggplot2"))
		suppressMessages(library("gridExtra"))
		suppressMessages(library("texreg"))
		suppressMessages(library("kableExtra"))
		suppressMessages(library("vtable"))
		suppressMessages(library("margins"))
		suppressMessages(library("purrr"))
		suppressMessages(library("interflex"))

		# Set Working Directory
		setwd('~/Dropbox/DHNR Philippine Police/3. EGAP (POP - One Sorsogon Reloaded)/3.14 - Replication Files/Effects on PNP/HRS_APSR_Replication/')

		# Define Functions
		'%!in%' <- function(x,y)!('%in%'(x,y))

	# =======================================================
	# Load Data
	# =======================================================	
		# Wide version of the officer survey (UoA: Officer)
		dt<-fread('2_Data/POLICE_wide.csv')
			# Create a dataframe version for some analyses
			df<-data.frame(dt)

		# Long version of the officer survey (UoA: Survey)
		dt_long<-fread('2_Data/POLICE_long.csv')	

		# Officer Embeddedness in the US
		oe<-fread('2_Data/USA_Embeddedness.csv')
			# Source: https://github.com/fivethirtyeight/data/tree/master/police-locals

		# Most Important Concerns
		issues<-fread('2_Data/MIP_issues.csv')


##################################################
############### 2. Main Analyses #################
##################################################

	# =======================================================
	# Average Treatment Effects - Panel Models
	# =======================================================	
		# --- Main Panel Models --- #
		outcomes<-c('top4match_pnplgu','officer_attitude_idx','trust_idx','empathy_idx','accountability_idx','corrupt_idx')
		estimates_list <- list()
		for (var in outcomes) {
		  estimates_list[[var]] <-
		    mod<-lm_robust(
		      as.formula(glue::glue("{ var }_e ~ Z_pop_assign + { var }_m + miss_bline")),
		      clusters = station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata",
		      data = dt) %>%
		    tidy 
		}
		estimates_panel <- rbindlist(estimates_list) 
		estimates_panel <- estimates_panel[term=='Z_pop_assign']
		estimates_panel$model<-'Panel'

		# --- Attach 90% confidence intervals --- # 
		ci<-matrix(nrow=6,ncol=2)
		m_top4_panel<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign + top4match_pnplgu_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[1,]<-confint(m_top4_panel,'Z_pop_assign',level=.9)[1:2]
		m_attitude_panel<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign + officer_attitude_idx_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[2,]<-confint(m_attitude_panel,'Z_pop_assign',level=.9)[1:2]
		m_trust_panel<-lm_robust(trust_idx_e ~ Z_pop_assign + trust_idx_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[3,]<-confint(m_trust_panel,'Z_pop_assign',level=.9)[1:2]
		m_empathy_panel<-lm_robust(empathy_idx_e ~ Z_pop_assign + empathy_idx_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[4,]<-confint(m_empathy_panel,'Z_pop_assign',level=.9)[1:2]
		m_account_panel<-lm_robust(accountability_idx_e ~ Z_pop_assign + accountability_idx_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[5,]<-confint(m_account_panel,'Z_pop_assign',level=.9)[1:2]
		m_corrupt_panel<-lm_robust(corrupt_idx_e ~ Z_pop_assign + corrupt_idx_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[6,]<-confint(m_corrupt_panel,'Z_pop_assign',level=.9)[1:2]
		cis<-data.table(ci90.low=ci[,1],ci90.high=ci[,2])
		estimates_panel<-cbind(estimates_panel,cis)

	# =======================================================
	# Average Treatment Effects - Cross-Sectional Models
	# =======================================================	
		# --- Main cross-sectional models --- # 
		estimates_list_cx <- list()
		for (var in outcomes) {
		  estimates_list_cx[[var]] <-
		    lm_robust(
		      as.formula(glue::glue("{ var }_e ~ Z_pop_assign ")),
		      clusters = station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata",
		      data = dt) %>%
		    tidy 
		}
		estimates_cross <- rbindlist(estimates_list_cx) 
		estimates_cross <- estimates_cross[term=='Z_pop_assign']
		estimates_cross$model<-'Endline Cross-Section'

		# --- Attach 90% confidence intervals --- # 
		ci<-matrix(nrow=6,ncol=2)
		m_top4_cross<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[1,]<-confint(m_top4_cross,'Z_pop_assign',level=.9)[1:2]
		m_attitude_cross<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[2,]<-confint(m_attitude_cross,'Z_pop_assign',level=.9)[1:2]
		m_trust_cross<-lm_robust(trust_idx_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[3,]<-confint(m_trust_cross,'Z_pop_assign',level=.9)[1:2]
		m_empathy_cross<-lm_robust(empathy_idx_e ~ Z_pop_assign,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[4,]<-confint(m_empathy_cross,'Z_pop_assign',level=.9)[1:2]
		m_account_cross<-lm_robust(accountability_idx_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[5,]<-confint(m_account_cross,'Z_pop_assign',level=.9)[1:2]
		m_corrupt_cross<-lm_robust(corrupt_idx_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[6,]<-confint(m_corrupt_cross,'Z_pop_assign',level=.9)[1:2]
		cis<-data.table(ci90.low=ci[,1],ci90.high=ci[,2])
		estimates_cross<-cbind(estimates_cross,cis)
		
		# --- Combine --- #
		estimates<-rbind(estimates_panel,estimates_cross)
		estimates[,model:=as.factor(model)]

	# =======================================================
	# Heterogeneous Effects - Embeddedness
	# =======================================================	
		# Regressions: Embeddedness -> Main Outcomes 
		m_top4_embed<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign*lives_in_sorsogon +top4match_pnplgu_m +miss_bline + rank_m_spo + station_m, data=df, clusters=station_m, se_type="stata")
		m_attitude_embed<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign*lives_in_sorsogon +officer_attitude_idx_m +miss_bline + rank_m_spo + station_m,data=df, clusters=station_m, se_type="stata")		
		m_trust_embed<-lm_robust(trust_idx_e ~ Z_pop_assign*lives_in_sorsogon +trust_idx_m +miss_bline + rank_m_spo + station_m, data=df, clusters=station_m, se_type="stata")
		m_empathy_embed<-lm_robust(empathy_idx_e ~ Z_pop_assign*lives_in_sorsogon +empathy_idx_m +miss_bline + rank_m_spo + station_m, data=df, clusters=station_m, se_type="stata")
		m_accountability_embed<-lm_robust(accountability_idx_e ~ Z_pop_assign*lives_in_sorsogon +accountability_idx_m +miss_bline + rank_m_spo + station_m, data=df, clusters=station_m, se_type="stata")
		m_corrupt_embed<-lm_robust(corrupt_idx_e ~ Z_pop_assign*lives_in_sorsogon +corrupt_idx_m +miss_bline + rank_m_spo + station_m, data=df, clusters=station_m, se_type="stata")

		# --- Margins: Main Outcomes --- #
		outcomes<-c('top4match_pnplgu','officer_attitude_idx','trust_idx','empathy_idx','accountability_idx','corrupt_idx')
		margins_embed_main <- list(); margins_embed_main90 <- list()
		for (var in outcomes) {
		  margins_embed_main[[var]] <- data.frame(summary(margins(lm_robust(as.formula(glue::glue("{ var }_e ~ Z_pop_assign + { var }_m + miss_bline + rank_m_spo + station_m + lives_in_sorsogon + Z_pop_assign*lives_in_sorsogon ")),clusters = station_m,se_type="stata",data = df,try_cholesky = TRUE),at=list(lives_in_sorsogon=0:1),variables='Z_pop_assign',vce="delta")),outcome=var)}

			# .90 CIs
			for (var in outcomes) {
		  	margins_embed_main90[[var]] <- data.frame(summary(margins(lm_robust(as.formula(glue::glue("{ var }_e ~ Z_pop_assign + { var }_m + miss_bline + rank_m_spo + station_m + lives_in_sorsogon + Z_pop_assign*lives_in_sorsogon ")),clusters = station_m,se_type="stata",data = df,try_cholesky = TRUE),at=list(lives_in_sorsogon=0:1),variables='Z_pop_assign',vce="delta"),level=.9),outcome=var)}

		 	# Clean
			margins_embed_main<-rbindlist(margins_embed_main)
			margins_embed_main90<-rbindlist(margins_embed_main90)
			margins_embed_main90<-margins_embed_main90[,c('lower','upper')]
			setnames(margins_embed_main90,c('lower','upper'),c('lower90','upper90'))
			margins_embed_main<-cbind(margins_embed_main,margins_embed_main90)
			margins_embed_main$label<-c('Share Citizen Concerns','Share Citizen Concerns','Officer Attitude Index','Officer Attitude Index','Trust','Trust','Empathy','Empathy','Accountability','Accountability','Corruption','Corruption')
			margins_embed_main$label <- factor (margins_embed_main$label,levels=c('Corruption','Accountability','Empathy','Trust','Officer Attitude Index','Share Citizen Concerns'))
			margins_embed_main$embedded<-factor(ifelse(margins_embed_main$lives_in_sorsogon==1,'Lives in Province','Lives Outside Province'))

	# =======================================================
	# Heterogeneous Effects - NPA (Rebel) Presence
	# =======================================================	

		# Regressions: NPA -> Main outcomes
			m_top4_npa<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +top4match_pnplgu_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")
			m_attitude_npa<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +officer_attitude_idx_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")
			m_trust_npa<-lm_robust(trust_idx_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_idx_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")		
			m_trust_alt_npa<-lm_robust(trust_idx_alt_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_idx_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")			
			m_empathy_npa<-lm_robust(empathy_idx_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +empathy_idx_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")
			m_account_npa<-lm_robust(accountability_idx_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +accountability_idx_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")
			m_corrupt_npa<-lm_robust(corrupt_idx_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +corrupt_idx_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")

		# Regressions: NPA -> Trust Sub-Indices & Sympathy
			m_tsafe_npa<-lm_robust(trust_safe_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_safe_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")
			m_tintent_npa<-lm_robust(trust_intent_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_intent_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")
			m_tinfo_npa<-lm_robust(trust_info_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_info_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")
			m_sympathy_npa<-lm_robust(npa_sympathy_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +npa_sympathy_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")

		# --- Margins: NPA --- # 
			# Knowledge
			margins_top4_npa<-data.frame(summary(margins(m_top4_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='top4match_pnplgu_e')
			margins_top4_npa90<-data.frame(summary(margins(m_top4_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='top4match_pnplgu_e')
			margins_top4_npa90<-margins_top4_npa90[,c('lower','upper')]
			setnames(margins_top4_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_top4_npa<-cbind(margins_top4_npa,margins_top4_npa90)

			# Attitude Index
			margins_attitude_npa<-data.frame(summary(margins(m_attitude_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='officer_attitude_idx_alt_e')
			margins_attitude_npa90<-data.frame(summary(margins(m_attitude_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='officer_attitude_idx_alt_e')
			margins_attitude_npa90<-margins_attitude_npa90[,c('lower','upper')]
			setnames(margins_attitude_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_attitude_npa<-cbind(margins_attitude_npa,margins_attitude_npa90)

			# Trust
			margins_trust_npa<-data.frame(summary(margins(m_trust_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='trust_idx_alt_e')
			margins_trust_npa90<-data.frame(summary(margins(m_trust_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='trust_idx_alt_e')
			margins_trust_npa90<-margins_trust_npa90[,c('lower','upper')]
			setnames(margins_trust_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_trust_npa<-cbind(margins_trust_npa,margins_trust_npa90)

			# Empathy
			margins_empathy_npa<-data.frame(summary(margins(m_empathy_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='empathy_idx_alt_e')
			margins_empathy_npa90<-data.frame(summary(margins(m_empathy_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='empathy_idx_alt_e')
			margins_empathy_npa90<-margins_empathy_npa90[,c('lower','upper')]
			setnames(margins_empathy_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_empathy_npa<-cbind(margins_empathy_npa,margins_empathy_npa90)

			# Accountability 
			margins_account_npa<-data.frame(summary(margins(m_account_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='accountability_idx_alt_e')
			margins_account_npa90<-data.frame(summary(margins(m_account_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='accountability_idx_alt_e')
			margins_account_npa90<-margins_account_npa90[,c('lower','upper')]
			setnames(margins_account_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_account_npa<-cbind(margins_account_npa,margins_account_npa90)

			# Corruption 
			margins_corrupt_npa<-data.frame(summary(margins(m_corrupt_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='corrupt_idx_alt_e')
			margins_corrupt_npa90<-data.frame(summary(margins(m_corrupt_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='corrupt_idx_alt_e')
			margins_corrupt_npa90<-margins_corrupt_npa90[,c('lower','upper')]
			setnames(margins_corrupt_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_corrupt_npa<-cbind(margins_corrupt_npa,margins_corrupt_npa90)

			# Trust - Safe
			margins_tsafe_npa<-data.frame(summary(margins(m_tsafe_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='trust_safe_e')
			margins_tsafe_npa90<-data.frame(summary(margins(m_tsafe_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='trust_safe_e')			
			margins_tsafe_npa90<-margins_tsafe_npa90[,c('lower','upper')]
			setnames(margins_tsafe_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_tsafe_npa<-cbind(margins_tsafe_npa,margins_tsafe_npa90)

			# Trust - Intent
			margins_tintent_npa<-data.frame(summary(margins(m_tintent_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='trust_intent_e')
			margins_tintent_npa90<-data.frame(summary(margins(m_tintent_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='trust_intent_e')
			margins_tintent_npa90<-margins_tintent_npa90[,c('lower','upper')]
			setnames(margins_tintent_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_tintent_npa<-cbind(margins_tintent_npa,margins_tintent_npa90)

			# Trust - Info
			margins_tinfo_npa<-data.frame(summary(margins(m_tinfo_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='trust_info_e')
			margins_tinfo_npa90<-data.frame(summary(margins(m_tinfo_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='trust_info_e')
			margins_tinfo_npa90<-margins_tinfo_npa90[,c('lower','upper')]
			setnames(margins_tinfo_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_tinfo_npa<-cbind(margins_tinfo_npa,margins_tinfo_npa90)

			# NPA Sympathy
			margins_sympathy_npa<-data.frame(summary(margins(m_sympathy_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='npa_sympathy_e')
			margins_sympathy_npa90<-data.frame(summary(margins(m_sympathy_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta"),level=.9),outcome='npa_sympathy_e')
			margins_sympathy_npa90<-margins_sympathy_npa90[,c('lower','upper')]
			setnames(margins_sympathy_npa90,c('lower','upper'),c('lower90','upper90'))
			margins_sympathy_npa<-cbind(margins_sympathy_npa,margins_sympathy_npa90)


#####################################################
############### 3. Robustness Tests #################
##################################################### 

	# =======================================================
	# Heterogeneous Effects - Embeddedness (With Controls)
	# =======================================================	
		# Individual controls + Station FEs  
		m_attitude_embed_ctrl1<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign*lives_in_sorsogon +officer_attitude_idx_m +miss_bline + station_m + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")
		 
		# Individual controls + Station Controls 
		m_attitude_embed_ctrl2<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign*lives_in_sorsogon +officer_attitude_idx_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin+ npa_2015_presence + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")
		
		# Individual controls (interacted) + Station Controls (interacted)
		m_attitude_embed_ctrl3<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign*lives_in_sorsogon +officer_attitude_idx_m +miss_bline + Z_pop_assign*c10_education + Z_pop_assign*c10_population + Z_pop_assign*c10_home_good +Z_pop_assign*c10_hhsize + Z_pop_assign*c10_catholic + Z_pop_assign*civ_crime_victim_idx_admin+ Z_pop_assign*npa_2015_presence + Z_pop_assign*rank_m_spo + Z_pop_assign*educ_high + Z_pop_assign*catholic + Z_pop_assign*age + Z_pop_assign*gender + Z_pop_assign*years_served_mun + Z_pop_assign*commute_m + Z_pop_assign*psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")

		# trust Index (with full battery of controls)
		m_trust_embed_ctrl3<-lm_robust(trust_idx_e ~ Z_pop_assign*lives_in_sorsogon +trust_idx_m +miss_bline + Z_pop_assign*c10_education + Z_pop_assign*c10_population + Z_pop_assign*c10_home_good +Z_pop_assign*c10_hhsize + Z_pop_assign*c10_catholic + Z_pop_assign*civ_crime_victim_idx_admin+ Z_pop_assign*npa_2015_presence + Z_pop_assign*rank_m_spo + Z_pop_assign*educ_high + Z_pop_assign*catholic + Z_pop_assign*age + Z_pop_assign*gender + Z_pop_assign*years_served_mun + Z_pop_assign*commute_m + Z_pop_assign*psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")

		# empathy Index (with full battery of controls)
		m_empathy_embed_ctrl3<-lm_robust(empathy_idx_e ~ Z_pop_assign*lives_in_sorsogon +empathy_idx_m +miss_bline + Z_pop_assign*c10_education + Z_pop_assign*c10_population + Z_pop_assign*c10_home_good +Z_pop_assign*c10_hhsize + Z_pop_assign*c10_catholic + Z_pop_assign*civ_crime_victim_idx_admin+ Z_pop_assign*npa_2015_presence + Z_pop_assign*rank_m_spo + Z_pop_assign*educ_high + Z_pop_assign*catholic + Z_pop_assign*age + Z_pop_assign*gender + Z_pop_assign*years_served_mun + Z_pop_assign*commute_m + Z_pop_assign*psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")

		# accountability Index (with full battery of controls)
		m_account_embed_ctrl3<-lm_robust(accountability_idx_e ~ Z_pop_assign*lives_in_sorsogon +accountability_idx_m +miss_bline + Z_pop_assign*c10_education + Z_pop_assign*c10_population + Z_pop_assign*c10_home_good +Z_pop_assign*c10_hhsize + Z_pop_assign*c10_catholic + Z_pop_assign*civ_crime_victim_idx_admin+ Z_pop_assign*npa_2015_presence + Z_pop_assign*rank_m_spo + Z_pop_assign*educ_high + Z_pop_assign*catholic + Z_pop_assign*age + Z_pop_assign*gender + Z_pop_assign*years_served_mun + Z_pop_assign*commute_m + Z_pop_assign*psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")

		# corrupt Index (with full battery of controls)
		m_corrupt_embed_ctrl3<-lm_robust(corrupt_idx_e ~ Z_pop_assign*lives_in_sorsogon +corrupt_idx_m +miss_bline + Z_pop_assign*c10_education + Z_pop_assign*c10_population + Z_pop_assign*c10_home_good +Z_pop_assign*c10_hhsize + Z_pop_assign*c10_catholic + Z_pop_assign*civ_crime_victim_idx_admin+ Z_pop_assign*npa_2015_presence + Z_pop_assign*rank_m_spo + Z_pop_assign*educ_high + Z_pop_assign*catholic + Z_pop_assign*age + Z_pop_assign*gender + Z_pop_assign*years_served_mun + Z_pop_assign*commute_m + Z_pop_assign*psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")

	# =======================================================
	# Heterogeneous Effects - NPA (With Controls)
	# =======================================================	
		# Station Controls
		m_trust_npa_ctrl1<-lm_robust(trust_idx_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_idx_m +miss_bline + rank_m_spo + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin,data=df,clusters=station_m, se_type="stata")
		 
		# Individual controls + Station Controls 
		m_trust_npa_ctrl2<-lm_robust(trust_idx_e ~ Z_pop_assign*npa_2015_presence +trust_idx_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin + rank_m_spo + lives_in_sorsogon + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")
			
		# Individual controls (interacted) + Station Controls (interacted)
		m_trust_npa_ctrl3<-lm_robust(trust_idx_e ~ Z_pop_assign*npa_2015_presence +trust_idx_m +miss_bline + Z_pop_assign*c10_education + Z_pop_assign*c10_population + Z_pop_assign*c10_home_good +Z_pop_assign*c10_hhsize + Z_pop_assign*c10_catholic + Z_pop_assign*civ_crime_victim_idx_admin + Z_pop_assign*rank_m_spo + Z_pop_assign*lives_in_sorsogon + Z_pop_assign*educ_high + Z_pop_assign*catholic + Z_pop_assign*age + Z_pop_assign*gender + Z_pop_assign*years_served_mun + Z_pop_assign*commute_m + Z_pop_assign*psm_idx_combined_m,data=df, clusters=station_m, se_type="stata")

		# Sub-index outcomes (with full battery of controls)
		m_tsafe_npa_ctrl2<-lm_robust(trust_safe_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_safe_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin + rank_m_spo + lives_in_sorsogon + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df,clusters=station_m, se_type="stata")
		m_tintent_npa_ctrl2<-lm_robust(trust_intent_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_intent_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin + rank_m_spo + lives_in_sorsogon + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df,clusters=station_m, se_type="stata")
		m_tinfo_npa_ctrl2<-lm_robust(trust_info_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_info_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin + rank_m_spo + lives_in_sorsogon + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df,clusters=station_m, se_type="stata")
		m_trust_alt_npa_ctrl2<-lm_robust(trust_idx_alt_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_idx_alt_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin + rank_m_spo + lives_in_sorsogon + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df,clusters=station_m, se_type="stata")		
		m_sympathy_npa_ctrl1<-lm_robust(npa_sympathy_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +npa_sympathy_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin + rank_m_spo,data=df,clusters=station_m, se_type="stata")
		m_sympathy_npa_ctrl2<-lm_robust(npa_sympathy_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +npa_sympathy_m +miss_bline + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + civ_crime_victim_idx_admin + rank_m_spo + lives_in_sorsogon + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df,clusters=station_m, se_type="stata")		

	# =======================================================
	# Panel Models - Alternate Outcome Index
	# =======================================================	
		# --- Main Panel Models --- #
		outcomes<-c('top4match_pnplgu','officer_attitude_idx_alt','trust_idx_alt','empathy_idx_alt','accountability_idx_alt','corrupt_idx_alt')
		estimates_list_alt <- list()
		for (var in outcomes) {
		  estimates_list_alt[[var]] <-
		    mod<-lm_robust(
		      as.formula(glue::glue("{ var }_e ~ Z_pop_assign + { var }_m + miss_bline")),
		      clusters = station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata",
		      data = dt) %>%
		    tidy 
		}
		estimates_panel_alt <- rbindlist(estimates_list_alt) 
		estimates_panel_alt <- estimates_panel_alt[term=='Z_pop_assign']
		estimates_panel_alt$model<-'Panel'

		# --- 90% confidence intervals --- # 
		ci<-matrix(nrow=6,ncol=2)
		m_top4_alt_panel<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign + top4match_pnplgu_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[1,]<-confint(m_top4_alt_panel,'Z_pop_assign',level=.9)[1:2]
		m_attitude_alt_panel<-lm_robust(officer_attitude_idx_alt_e ~ Z_pop_assign + officer_attitude_idx_alt_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[2,]<-confint(m_attitude_panel,'Z_pop_assign',level=.9)[1:2]
		m_trust_alt_panel<-lm_robust(trust_idx_alt_e ~ Z_pop_assign + trust_idx_alt_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[3,]<-confint(m_trust_alt_panel,'Z_pop_assign',level=.9)[1:2]
		m_empathy_alt_panel<-lm_robust(empathy_idx_alt_e ~ Z_pop_assign + empathy_idx_alt_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[4,]<-confint(m_empathy_alt_panel,'Z_pop_assign',level=.9)[1:2]
		m_account_alt_panel<-lm_robust(accountability_idx_alt_e ~ Z_pop_assign + accountability_idx_alt_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[5,]<-confint(m_account_alt_panel,'Z_pop_assign',level=.9)[1:2]
		m_corrupt_alt_panel<-lm_robust(corrupt_idx_alt_e ~ Z_pop_assign + corrupt_idx_alt_m +miss_bline,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[6,]<-confint(m_corrupt_alt_panel,'Z_pop_assign',level=.9)[1:2]
		cis<-data.table(ci90.low=ci[,1],ci90.high=ci[,2])
		estimates_panel_alt<-cbind(estimates_panel_alt,cis)

	# =======================================================
	# Cross-Sectional Models - Alternate Outcome Index
	# =======================================================	
		# --- Main cross-sectional models --- # 
		estimates_list_cx_alt <- list()
		for (var in outcomes) {
		  estimates_list_cx_alt[[var]] <-
		    lm_robust(
		      as.formula(glue::glue("{ var }_e ~ Z_pop_assign ")),
		      clusters = station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata",
		      data = dt) %>%
		    tidy
		}
		estimates_cross_alt <- rbindlist(estimates_list_cx_alt) 
		estimates_cross_alt <- estimates_cross_alt[term=='Z_pop_assign']
		estimates_cross_alt$model<-'Endline Cross-Section'

		# --- 90% confidence intervals --- # 
		ci<-matrix(nrow=6,ncol=2)
		m_top4_alt_cross<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[1,]<-confint(m_top4_alt_cross,'Z_pop_assign',level=.9)[1:2]
		m_attitude_alt_cross<-lm_robust(officer_attitude_idx_alt_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[2,]<-confint(m_attitude_alt_cross,'Z_pop_assign',level=.9)[1:2]
		m_trust_alt_cross<-lm_robust(trust_idx_alt_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[3,]<-confint(m_trust_alt_cross,'Z_pop_assign',level=.9)[1:2]
		m_empathy_alt_cross<-lm_robust(empathy_idx_alt_e ~ Z_pop_assign,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[4,]<-confint(m_empathy_alt_cross,'Z_pop_assign',level=.9)[1:2]
		m_account_alt_cross<-lm_robust(accountability_idx_alt_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[5,]<-confint(m_account_alt_cross,'Z_pop_assign',level=.9)[1:2]
		m_corrupt_alt_cross<-lm_robust(corrupt_idx_alt_e ~ Z_pop_assign ,data=dt,clusters=station_m, fixed_effects=~station_m+rank_m_spo,se_type="stata"); ci[6,]<-confint(m_corrupt_alt_cross,'Z_pop_assign',level=.9)[1:2]
		cis<-data.table(ci90.low=ci[,1],ci90.high=ci[,2])
		estimates_cross_alt<-cbind(estimates_cross_alt,cis)
		
		# --- Combine --- #
		estimates_alt<-rbind(estimates_panel_alt,estimates_cross_alt)
		estimates_alt[,model:=as.factor(model)]

	# =======================================================
	# Heterogeneous Effects - Alternate Outcome Index
	# =======================================================	
		# --- Embeddedness --- #
			# Regressions: Embeddedness -> Officer Attitude Index 
			m_attitude_alt_embed<-lm_robust(officer_attitude_idx_alt_e ~ Z_pop_assign*lives_in_sorsogon +officer_attitude_idx_alt_m +miss_bline + rank_m_spo + station_m,data=df, clusters=station_m, se_type="stata")
			
			# Margins: Officer Attitude Index
			margins_embed_alt_main<-data.frame(summary(margins(m_attitude_alt_embed,at=list(lives_in_sorsogon=0:1),variables='Z_pop_assign',vce="delta")),outcome='officer_attitude_idx_alt')
			margins_embed_alt_main$label<-factor(c('Officer Attitude Index (alt)','Officer Attitude Index (alt)'))
			margins_embed_alt_main$embedded<-factor(ifelse(margins_embed_alt_main$lives_in_sorsogon==1,'Lives in Province','Lives Outside Province'))

		# --- NPA --- #
			# Regressions: NPA -> Trust & Sympathy
			m_trustidx_alt_npa<-lm_robust(trust_idx_alt_e ~ Z_pop_assign + npa_2015_presence + Z_pop_assign*npa_2015_presence +trust_idx_alt_m +miss_bline + rank_m_spo,data=df,clusters=station_m, se_type="stata")

			# Margins: NPA -> Trust & Sympathy
			margins_trustidx_alt_npa<-data.frame(summary(margins(m_trustidx_alt_npa,at=list(npa_2015_presence=as.numeric(names(table(dt$npa_2015_presence)))),variables='Z_pop_assign',vce="delta")),outcome='trust_idx_alt_e')	

	# =======================================================
	# Spillover Analysis 
	# =======================================================

		# --- Het effects: station size --- # 
			m_attitude_statsize<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign*station_size + rank_m_spo + miss_bline + officer_attitude_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic,data=df, clusters=station_m, se_type="stata")
			m_top4_statsize<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign*station_size + rank_m_spo + miss_bline + top4match_pnplgu_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic,data=df, clusters=station_m, se_type="stata")
			m_trust_statsize<-lm_robust(trust_idx_e ~ Z_pop_assign*station_size + rank_m_spo + miss_bline + trust_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic,data=df, clusters=station_m, se_type="stata")
			m_empathy_statsize<-lm_robust(empathy_idx_e ~ Z_pop_assign*station_size + rank_m_spo + miss_bline + empathy_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic,data=df, clusters=station_m, se_type="stata")
			m_accountability_statsize<-lm_robust(accountability_idx_e ~ Z_pop_assign*station_size + rank_m_spo + miss_bline + accountability_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic,data=df, clusters=station_m, se_type="stata")
			m_corruption_statsize<-lm_robust(corrupt_idx_e ~ Z_pop_assign*station_size + rank_m_spo + miss_bline + corrupt_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic,data=df, clusters=station_m, se_type="stata")

		# --- Het effects: station density --- # 
			m_attitude_statdensity<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign*density_station + rank_m_spo + miss_bline + officer_attitude_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + station_size,data=df, clusters=station_m, se_type="stata")
			m_top4_statdensity<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign*density_station + rank_m_spo + miss_bline + top4match_pnplgu_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + station_size,data=df, clusters=station_m, se_type="stata")
			m_trust_statdensity<-lm_robust(trust_idx_e ~ Z_pop_assign*density_station + rank_m_spo + miss_bline + trust_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + station_size,data=df, clusters=station_m, se_type="stata")
			m_empathy_statdensity<-lm_robust(empathy_idx_e ~ Z_pop_assign*density_station + rank_m_spo + miss_bline + empathy_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + station_size,data=df, clusters=station_m, se_type="stata")
			m_accountability_statdensity<-lm_robust(accountability_idx_e ~ Z_pop_assign*density_station + rank_m_spo + miss_bline + accountability_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + station_size,data=df, clusters=station_m, se_type="stata")
			m_corruption_statdensity<-lm_robust(corrupt_idx_e ~ Z_pop_assign*density_station + rank_m_spo + miss_bline + corrupt_idx_m + npa_2015_presence + civ_crime_victim_idx_admin + c10_education + c10_population + c10_home_good +c10_hhsize + c10_catholic + station_size,data=df, clusters=station_m, se_type="stata")

		# --- Het effects: Direct Ties to Treatment --- # 
			# Knowledge
			m_knowledge_spillany<-lm_robust(top4match_pnplgu_e ~ spill_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0,], clusters=station_e, se_type="stata")
			m_knowledge_spillfam<-lm_robust(top4match_pnplgu_e ~ spill_family_direct + rank_m_spo + station_e,data=df[df$Z_pop_assign==0,], clusters=station_e, se_type="stata")
			m_knowledge_spillbeat<-lm_robust(top4match_pnplgu_e ~ spill_beat_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0,], clusters=station_e, se_type="stata")

			# Attitude Index
			m_attitude_spillany<-lm_robust(officer_attitude_idx_e ~ spill_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0 & df$station_e!='',], clusters=station_e, se_type="stata")
			m_attitude_spillfam<-lm_robust(officer_attitude_idx_e ~ spill_family_direct + rank_m_spo + station_e,data=df[df$Z_pop_assign==0 & df$station_e!='',], clusters=station_e, se_type="stata")
			m_attitude_spillbeat<-lm_robust(officer_attitude_idx_e ~ spill_beat_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0 & df$station_e!='',], clusters=station_e, se_type="stata")

			# Sub-Indices
			m_trust_spillany<-lm_robust(trust_idx_e ~ spill_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0 & df$station_e!='',], clusters=station_e, se_type="stata")
			m_empathy_spillany<-lm_robust(empathy_idx_e ~ spill_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0 & df$station_e!='',], clusters=station_e, se_type="stata")
			m_accountability_spillany<-lm_robust(accountability_idx_e ~ spill_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0 & df$station_e!='',], clusters=station_e, se_type="stata")
			m_corrupt_spillany<-lm_robust(corrupt_idx_e ~ spill_any + rank_m_spo + station_e,data=df[df$Z_pop_assign==0 & df$station_e!='',], clusters=station_e, se_type="stata")

	# =======================================================
	# Alternative "Shared Concerns" Measures
	# =======================================================
		# Using top 3 match instead of top 4 
		m_top3_panel<-lm_robust(top3match_pnplgu_e ~ Z_pop_assign + top3match_pnplgu_m + miss_bline,clusters=station_m,fixed_effect=~station_m+rank_m_spo,se_type='stata',data=dt)
		m_top3_cross<-lm_robust(top3match_pnplgu_e ~ Z_pop_assign ,clusters=station_m,fixed_effect=~station_m+rank_m_spo,se_type='stata',data=dt)

		# Using survey question on what officers thought citizens thought
		m_top4know_panel<-lm_robust(top4match_civlgu_e ~ Z_pop_assign + top4match_civlgu_m + miss_bline,clusters=station_m,fixed_effect=~station_m+rank_m_spo,se_type='stata',data=dt)
		m_top4know_cross<-lm_robust(top4match_civlgu_e ~ Z_pop_assign ,clusters=station_m,fixed_effect=~station_m+rank_m_spo,se_type='stata',data=dt)

		# Using match with ordinary citizens 
		m_top4civ_panel<-lm_robust(top4match_pnpciv_e ~ Z_pop_assign + top4match_pnpciv_m + miss_bline,clusters=station_m,fixed_effect=~station_m+rank_m_spo,se_type='stata',data=dt)
		m_top4civ_cross<-lm_robust(top4match_pnpciv_e ~ Z_pop_assign ,clusters=station_m,fixed_effect=~station_m+rank_m_spo,se_type='stata',data=dt)			

	# =======================================================
	# Heterogeneous Effects - Alternate Embeddedness Measure
	# =======================================================	
		# Lives in Municipality 
		m_top4_embedmun<-lm_robust(top4match_pnplgu_e ~ Z_pop_assign*lives_where_works +top4match_pnplgu_m +miss_bline + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, fixed_effects = ~ station_m, clusters=station_m, se_type="stata")			
		m_attitude_embedmun<-lm_robust(officer_attitude_idx_e ~ Z_pop_assign*lives_where_works +officer_attitude_idx_m +miss_bline + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, fixed_effects = ~ station_m, clusters=station_m, se_type="stata")
		m_trust_embedmun<-lm_robust(trust_idx_e ~ Z_pop_assign*lives_where_works +trust_idx_m +miss_bline + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, fixed_effects = ~ station_m, clusters=station_m, se_type="stata")
		m_empathy_embedmun<-lm_robust(empathy_idx_e ~ Z_pop_assign*lives_where_works +empathy_idx_m +miss_bline + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, fixed_effects = ~ station_m, clusters=station_m, se_type="stata")
		m_accountability_embedmun<-lm_robust(accountability_idx_e ~ Z_pop_assign*lives_where_works +accountability_idx_m +miss_bline + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, fixed_effects = ~ station_m, clusters=station_m, se_type="stata")
		m_corrupt_embedmun<-lm_robust(corrupt_idx_e ~ Z_pop_assign*lives_where_works +corrupt_idx_m +miss_bline + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, fixed_effects = ~ station_m, clusters=station_m, se_type="stata")
		m_attitude_alt_embedmun<-lm_robust(officer_attitude_idx_alt_e ~ Z_pop_assign*lives_where_works +officer_attitude_idx_alt_m +miss_bline + rank_m_spo + educ_high + catholic + age + gender + years_served_mun + commute_m + psm_idx_combined_m,data=df, fixed_effects = ~ station_m, clusters=station_m, se_type="stata")	
	
	
#########################################################
############### 4. Figures - Main Paper #################
######################################################### 

	# =======================================================
	# Main outcomes
	# =======================================================
		# Prepare dataset for figure
		gg_estimates<-as.data.table(estimates)
		gg_estimates[,p_print:=round(p.value,3)]
		vars<-c('Share Local Concerns','Officer Attitude Index','Trust','Empathy','Accountability','Corruption')
		gg_estimates$label<-c(vars,vars)
		gg_estimates$label <- factor(gg_estimates$label, levels = c('Corruption','Accountability','Empathy','Trust','Officer Attitude Index','Share Local Concerns'))

		# Figure
		pdf('4_Figures/Fig1_Main_CoefficientPlot.pdf',width=8,height=6)
		ggplot(gg_estimates,aes(x=estimate, y=label, shape = model, color = model, fill = model)) + 
		geom_point(size=3,position = position_dodge(width=.35)) + 
	  	geom_errorbarh(aes(xmin=conf.low, xmax=conf.high), position=position_dodge(width=.35), height=0) + 
	  	geom_errorbarh(aes(xmin=ci90.low, xmax=ci90.high), position=position_dodge(width=.35), height=0,size=1.2) + 
	  	geom_vline(xintercept=0, linetype="dotted") +
		scale_color_manual(values = c("#D7C49EFF","#343148FF"), name = "") + 
		annotate('rect', xmin=-.36, xmax=.45, ymin=0, ymax=4.4, alpha=.1, fill='cornflowerblue') +
	  	scale_fill_manual(values = c("#D7C49EFF","#343148FF"), name = "") + 
	  	scale_shape_manual(values = c(19, 23)) + 
		labs(title = "", x = "Estimated treatment effect", y = "", colour = "", shape = "") +
		theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90"), legend.position='bottom', axis.title=element_text(face = "bold",size=14),legend.text=element_text(size=12),legend.title=element_text(face='bold',size=13),axis.text.x=element_text(size=12),axis.text.y=element_text(size=13)) 
		dev.off()

	# =======================================================
	# Marginal Effects - Embeddedness
	# =======================================================	
		# Figure
		pdf('4_Figures/Fig2_Margins_Embedded_Main.pdf',width=8,height=6)
		ggplot(margins_embed_main,aes(x=AME, y=label, shape = embedded, color = embedded, fill = embedded)) + 
		xlim(-.52,.86) +
		geom_point(size=3,position = position_dodge(width=.35)) + 
	  	geom_errorbarh(aes(xmin=lower, xmax=upper), position=position_dodge(width=.35), height=0) + 
  		geom_errorbarh(aes(xmin=lower90, xmax=upper90), position=position_dodge(width=.35), height=0,size=1.2) + 
	  	geom_vline(xintercept=0, linetype="dotted") +
			scale_color_manual(values = c("#D7C49EFF","#191970"), name = "") + 
	  	scale_fill_manual(values = c("#D7C49EFF","#191970"), name = "") + 
	  	scale_shape_manual(values = c(19, 23)) + 
		annotate('rect', xmin=-.52, xmax=.86, ymin=0, ymax=4.4, alpha=.1, fill='cornflowerblue') +
		labs(title = "", x = "Estimated treatment effect", y = "", colour = "", shape = "") +
		theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90"), legend.position='bottom', axis.title=element_text(face = "bold",size=14),legend.text=element_text(size=12),legend.title=element_text(face='bold',size=13),axis.text.x=element_text(size=12),axis.text.y=element_text(size=13)) 
		dev.off()


	# =======================================================
	# US Officer Embeddedness
	# =======================================================
		# Prep
		oe<-oe[order(-Percentage)]
		oe$City<-factor(oe$City,levels=rev(oe$City))

		# Export Plot
		pdf('4_Figures/Fig3_USEmbeddedness.pdf',width=10,height=7)
		ggplot(oe, aes(x=Percentage, y=City,size=Size)) + geom_point() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90"), legend.position='bottom', axis.title=element_text(face = "bold",size=14),legend.text=element_text(size=12),legend.title=element_text(face='bold',size=13),axis.text.x=element_text(size=12),axis.text.y=element_text(size=13)) + xlab("Percentage of Officers Who Live in the City") + ylab("") + guides(alpha = "none", size = "none") +	xlim(0,100) 
		dev.off()

	# =======================================================
	# Marginal Effects - NPA
	# =======================================================		
		# --- NPA -> Trust Index --- # 
			p_trust_npa<-ggplot(margins_trust_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + ylim(-.75,.75) +
			ggtitle("Trust Index") + xlab("") + ylab("Treatment Effect") +
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(6,15,15,6),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- NPA -> Trust Safety --- # 
			p_tsafe_npa<-ggplot(margins_tsafe_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + ylim(-30,15) +		
			ggtitle("Trust (Safety)") + xlab("") + ylab("Treatment Effect") +
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(6,12,6,6),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- NPA -> Trust Intent --- # 
			p_tintent_npa<-ggplot(margins_tintent_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + ylim(-30,15) +	
			ggtitle("Trust (Intent)") + xlab("NPA Presence") + ylab("") +							
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(6,12,6,6),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- NPA -> Trust Information --- # 
			p_tinfo_npa<-ggplot(margins_tinfo_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + ylim(-30,15) +
			ggtitle("Trust (Information)") +	xlab("") + ylab("") +
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(6,12,6,6),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- NPA -> NPA Sympathy --- # 
			p_sympathy_npa<-ggplot(margins_sympathy_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + ylim(-25,25) +
			ggtitle("Perceived NPA Sympathy") +	xlab("") + ylab("") +			
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(6,12,6,6),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- Export Combined NPA Marginal Effects Plot: TRUST --- # 
			pdf('4_Figures/Fig4_Margins_NPA_trust.pdf',width=13,height=11)
			grid.arrange(p_trust_npa, p_sympathy_npa, p_tsafe_npa, p_tintent_npa, p_tinfo_npa, nrow = 2, heights=c(1.5,1),
             layout_matrix = rbind(c(1,1,1,2,2,2), c(3,3,4,4,5,5)))
			dev.off()


#######################################################
############### 5. Figures - Appendix #################
####################################################### 

	# =======================================================
	# Descriptive - Outcome Variables Distribution
	# =======================================================
		# --- Attitude Indices --- #
			# Prep
			y<-dt[!is.na(sid_m),c('officer_attitude_idx_raw_m','trust_idx_raw_m','empathy_idx_raw_m','accountability_idx_raw_m','corrupt_idx_raw_m')]
			plotdata_raw<- data.table(melt(y))
			plotdata_raw[variable=='officer_attitude_idx_raw_m',variable:='Officer Attitude Index']
			plotdata_raw[variable=='trust_idx_raw_m',variable:='Trust Index']
			plotdata_raw[variable=='empathy_idx_raw_m',variable:='Empathy Index']
			plotdata_raw[variable=='accountability_idx_raw_m',variable:='Accountability Index']
			plotdata_raw[variable=='corrupt_idx_raw_m',variable:='Corruption Index']
			setnames(plotdata_raw,'variable','Variable')  

			# Plot 
			distribution_attitudes<-ggplot(plotdata_raw, aes(x=value, fill=Variable)) + xlim(0,100) + 
			theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90"), legend.position=c(0.3,0.8), axis.title=element_text(face = "bold",size=14),legend.text=element_text(size=14),legend.title=element_text(face='bold',size=15),axis.text.x=element_text(size=12),axis.text.y=element_text(size=12)) + labs(x='Attitude Index Values',y='Density') + 
			geom_density(aes(group=Variable, colour=Variable, fill=Variable), alpha=0.2, bw=6) + scale_color_brewer(palette = "Set1") + scale_fill_brewer(palette = "Set1")  

		# --- Shared Concerns --- # 
			tab<-as.data.table(table(dt[!is.na(sid_m)]$top4match_pnplgu_m))
			distribution_knowledge<-ggplot(data=tab, aes(x=V1, y=N)) +
			geom_bar(stat="identity", fill="steelblue")+
			theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90"),axis.title=element_text(face = "bold",size=14),axis.text.x=element_text(size=12),axis.text.y=element_text(size=12)) + labs(x='Shared Concerns',y='Frequency') 

		# --- Export Combined Plot --- # 
			pdf('4_Figures/FigA1_Distribution_Outcomes.pdf',width=11,height=5.5)
			grid.arrange(distribution_knowledge, distribution_attitudes, nrow = 1,widths=c(.58,1))
			dev.off()

	# =======================================================
	# Descriptive - Treatment Compliance/Intensity
	# =======================================================	
		# Meetings Attended
		pdf('4_Figures/FigA2.1_Histogram_Attend1.pdf',width=7,height=3.5) 
	    hist(dt[!is.na(sid_e)]$attended_total,main="Community Policing Meetings Attended",xlab="Meetings Attended", xlim=c(0,10),freq=TRUE,breaks=10) 
	    dev.off()

	    # Attendance Rate
		pdf('4_Figures/FigA2.2_Histogram_Attend2.pdf',width=7,height=3.5) 
	    hist(dt[!is.na(sid_e)]$attend_rate,main="Community Policing Attendance Rate",xlab="Attendance Rate", xlim=c(0,1),freq=FALSE,breaks=10) 
	    dev.off()

	# =======================================================
	# Marginal Effects - NPA (Main Outcomes)
	# =======================================================	
		# --- NPA -> Knowledge --- # 
			p_top4_npa<-ggplot(margins_top4_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + ylim(-1,1) +
			ggtitle("Share Citizen Concerns") + xlab("NPA Presence") + ylab("Treatment Effect") +
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(6,15,15,6),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- NPA -> Attitude Index --- # 
			p_attitude_npa<-ggplot(margins_attitude_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + ylim(-1,1) +
			ggtitle("Officer Attitude Index") + xlab("NPA Presence") + ylab("") +
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(6,15,15,6),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- NPA -> Trust Index --- # 
			p_trust_npa2<-ggplot(margins_trust_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + 
			ggtitle("Trust") + xlab("") + ylab("") + 
			scale_x_continuous(n.breaks=3) + scale_y_continuous(n.breaks=4,limits=c(-.75,.5)) +
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(2,2,0,-25),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3)))

		# --- NPA -> Empathy Index --- # 
			p_empathy_npa<-ggplot(margins_empathy_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + 
			ggtitle("Empathy") + xlab("NPA Presence") + ylab("") +
			scale_x_continuous(n.breaks=3) + scale_y_continuous(n.breaks=4,limits=c(-.75,.5)) +			
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(2,2,0,-25),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3))) 	

		# --- NPA -> Accountability Index --- # 
			p_account_npa<-ggplot(margins_account_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + 
			ggtitle("Accountability") + xlab("NPA Presence") + ylab("") +
			scale_x_continuous(n.breaks=3) + scale_y_continuous(n.breaks=4,limits=c(-.75,.5)) +			
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(2,2,0,-25),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3))) 	

		# --- NPA -> Corruption Index --- # 
			p_corrupt_npa<-ggplot(margins_corrupt_npa, aes(x = npa_2015_presence)) + 
			geom_line(aes(y = AME)) +
			geom_ribbon(aes(ymin = lower, ymax = upper), alpha = .12) + 
			geom_ribbon(aes(ymin = lower90, ymax = upper90), alpha = 0.2) + 
			geom_hline(yintercept = 0, linetype = "dashed") + 
			ggtitle("Corruption") + xlab("") + ylab("") +			
			scale_x_continuous(n.breaks=3) + scale_y_continuous(n.breaks=4,limits=c(-.75,.5)) +
			theme(plot.title = element_text(hjust = 0.5, size=20, face='bold')) +       
	        theme(panel.grid.major = element_line(colour="grey90", size=0.1), panel.grid.minor = element_blank(), panel.background = element_blank(), axis.line = element_line(colour = "grey90")) + 
	        theme(axis.title=element_text(size=24),legend.text=element_text(size=24),legend.title=element_text(face='bold',size=16),axis.text.x=element_text(size=16),axis.text.y=element_text(size=20)) + 
	        theme(axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) + theme(plot.title = element_text(margin = margin(t = 0, b = 45))) +
	        theme(plot.margin=unit(c(2,4,0,-25),"points")) +
	      	guides(colour = guide_legend(override.aes = list(size=3))) 	

		# --- Export Combined NPA Marginal Effects Plot: MAIN OUTCOMES --- # 
			pdf('4_Figures/FigA3_Margins_NPA_main.pdf',width=13,height=11)
			grid.arrange(p_top4_npa, p_attitude_npa, p_trust_npa2, p_empathy_npa, p_account_npa,p_corrupt_npa, nrow = 2, heights=c(1.5,1),
             layout_matrix = rbind(c(1,1,1,1,2,2,2,2), c(3,3,4,4,5,5,6,6)))
			dev.off()

	# =======================================================
	# Marginal Effects - NPA (Interflex)
	# =======================================================	
		# Outcome: Trust Index
		p_trustidx_npa_flex <- interflex(Y = "trust_idx_alt_e", D = "Z_pop_assign", X = "npa_2015_presence", Z = c("trust_idx_alt_m","miss_bline","rank_m_spo"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Trust Index",xlab='',ylab='Treatment Effect',xlim=c(0,1.25),ylim=c(-.75,.75),theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=5)

		# Outcome: Trust Safe
		p_tsafe_npa_flex <- interflex(Y = "trust_safe_e", D = "Z_pop_assign", X = "npa_2015_presence", Z = c("trust_safe_m","miss_bline","rank_m_spo"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Trust (Safety)",xlab='',ylab='Treatment Effect',xlim=c(0,1.25),ylim=c(-30,15),theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=5)

		# Outcome: Trust Intent
		p_tintent_npa_flex <- interflex(Y = "trust_intent_e", D = "Z_pop_assign", X = "npa_2015_presence", Z = c("trust_intent_m","miss_bline","rank_m_spo"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Trust (Intent)",xlab='NPA Presence',ylab='',xlim=c(0,1.25),ylim=c(-30,15),theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=5)

		# Outcome: Trust Info
		p_tinfo_npa_flex <- interflex(Y = "trust_info_e", D = "Z_pop_assign", X = "npa_2015_presence", Z = c("trust_info_m","miss_bline","rank_m_spo"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Trust (Information)",xlab='',ylab='',xlim=c(0,1.25),ylim=c(-30,15),theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=5)

		# Outcome: NPA Sympathy
		p_sympathy_npa_flex <- interflex(Y = "npa_sympathy_e", D = "Z_pop_assign", X = "npa_2015_presence", Z = c("npa_sympathy_m","miss_bline","rank_m_spo"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Perceived NPA Sympathy",xlab='',ylab='',xlim=c(0,1.25),ylim=c(-25,25),theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=5)		

		# --- Export Combined NPA Marginal Effects Plot: TRUST SUB-INDICES --- # 
		pdf('4_Figures/FigA4_Margins_Interflex_NPA.pdf',width=13,height=11)
		grid.arrange(p_trustidx_npa_flex$figure, p_sympathy_npa_flex$figure, p_tsafe_npa_flex$figure, p_tintent_npa_flex$figure, p_tinfo_npa_flex$figure, nrow = 2, heights=c(1.5,1),
         layout_matrix = rbind(c(1,1,1,2,2,2), c(3,3,4,4,5,5)))
		dev.off()			

	# =======================================================
	# Spillovers Marginal Effects (Interflex)
	# =======================================================
		# --- Station Network Density --- # 
	      p_top4_statdensity <- interflex(Y = "top4match_pnplgu_e", D = "Z_pop_assign", X = "density_station", Z = c("rank_m_spo","miss_bline","top4match_pnplgu_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic","station_size"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Share Concerns",xlab='Station Network Density',ylab='Treatment Effect',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_attitude_statdensity <- interflex(Y = "officer_attitude_idx_e", D = "Z_pop_assign", X = "density_station", Z = c("rank_m_spo","miss_bline","officer_attitude_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic","station_size"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Attitude Index",xlab='Station Network Density',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_trust_statdensity <- interflex(Y = "trust_idx_e", D = "Z_pop_assign", X = "density_station", Z = c("rank_m_spo","miss_bline","trust_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic","station_size"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Trust",xlab='',ylab='Treatment Effect',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_empathy_statdensity <- interflex(Y = "empathy_idx_e", D = "Z_pop_assign", X = "density_station", Z = c("rank_m_spo","miss_bline","empathy_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic","station_size"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Empathy",xlab='Station Density',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_account_statdensity <- interflex(Y = "accountability_idx_e", D = "Z_pop_assign", X = "density_station", Z = c("rank_m_spo","miss_bline","accountability_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic","station_size"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Accountability",xlab='Station Density',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_corrupt_statdensity <- interflex(Y = "corrupt_idx_e", D = "Z_pop_assign", X = "density_station", Z = c("rank_m_spo","miss_bline","corrupt_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic","station_size"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Corruption",xlab='',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)	     

	   # --- Export Plots --- # 
			pdf('4_Figures/FigA5_Margins_Spillover_Density.pdf',width=13,height=11)
			grid.arrange(p_top4_statdensity$figure, p_attitude_statdensity$figure, p_trust_statdensity$figure, p_empathy_statdensity$figure, p_account_statdensity$figure, p_corrupt_statdensity$figure, nrow = 2, heights=c(1.5,1),layout_matrix = rbind(c(1,1,1,1,2,2,2,2), c(3,3,4,4,5,5,6,6)))
			dev.off()

		# --- Station Size --- # 
	      p_top4_statsize <- interflex(Y = "top4match_pnplgu_e", D = "Z_pop_assign", X = "station_size", Z = c("rank_m_spo","miss_bline","top4match_pnplgu_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Share Concerns",xlab='Station Size',ylab='Treatment Effect',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_attitude_statsize <- interflex(Y = "officer_attitude_idx_e", D = "Z_pop_assign", X = "station_size", Z = c("rank_m_spo","miss_bline","officer_attitude_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Attitude Index",xlab='Station Size',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_trust_statsize <- interflex(Y = "trust_idx_e", D = "Z_pop_assign", X = "station_size", Z = c("rank_m_spo","miss_bline","trust_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Trust",xlab='',ylab='Treatment Effect',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_empathy_statsize <- interflex(Y = "empathy_idx_e", D = "Z_pop_assign", X = "station_size", Z = c("rank_m_spo","miss_bline","empathy_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Empathy",xlab='Station Size',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_account_statsize <- interflex(Y = "accountability_idx_e", D = "Z_pop_assign", X = "station_size", Z = c("rank_m_spo","miss_bline","accountability_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Accountability",xlab='Station Size',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)
	      p_corrupt_statsize <- interflex(Y = "corrupt_idx_e", D = "Z_pop_assign", X = "station_size", Z = c("rank_m_spo","miss_bline","corrupt_idx_m","npa_2015_presence","civ_crime_victim_idx_admin","c10_education","c10_population","c10_home_good","c10_hhsize","c10_catholic"), data = df, estimator = "binning", vcov.type = "cluster", cl="station_m",bin.labs = FALSE,Xdistr = "density",main = "Corruption",xlab='',ylab='',theme.bw = TRUE,cex.main=1.4,cex.lab=1.4,cex.axis=1.2,na.rm=TRUE,nbins=3)	     

	   # --- Export Plots --- # 
			pdf('4_Figures/FigA6_Margins_Spillover_Size.pdf',width=13,height=11)
			grid.arrange(p_top4_statsize$figure, p_attitude_statsize$figure, p_trust_statsize$figure, p_empathy_statsize$figure, p_account_statsize$figure, p_corrupt_statsize$figure, nrow = 2, heights=c(1.5,1),layout_matrix = rbind(c(1,1,1,1,2,2,2,2), c(3,3,4,4,5,5,6,6)))
			dev.off()


############################################
############### 6. Tables  #################
############################################

	# =======================================================
	# Descriptive - Most important concerns
	# =======================================================
		# Prep
		issues_table<-data.table(Issue=c('Public Intoxication','Theft','Illegal Gambling','Police Abuse','Sexual Harassment','Robbery','Vehicle Theft','Vehicle Accidents','Illegal Guns','Domestic Abuse','Murder','Illegal Drug Use','Rape'),PNP=c(round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_intox)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_theft)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_gambling)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_polabuse)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_sexharass)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_robbery)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_carnap)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_vehicle)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_firearms)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_domabuse)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_murder)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_drugs)),1),round(100*mean(na.omit(issues[respondent_cat=='officer']$top3_rape)),1)),Citizens=c(round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_intox)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_theft)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_gambling)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_polabuse)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_sexharass)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_robbery)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_carnap)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_vehicle)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_firearms)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_domabuse)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_murder)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_drugs)),1),round(100*mean(na.omit(issues[respondent_cat=='citizen']$top3_rape)),1)))
		issues_table[,Difference:=PNP-Citizens]
		issues_table[,PNP:=paste0(PNP,"\\%")]; issues_table[,Citizens:=paste0(Citizens,"\\%")]; issues_table[,Difference:=paste0(Difference,"\\%")]

      # Export Table
		kable(issues_table, digits =1,align = "llll",format='latex',position='h') %>%
		save_kable('3_Tables/0_raw/Tab1_SummaryConcerns_raw.tex')

	# =======================================================
	# Descriptive - Summary Statistics Table
	# =======================================================
		# Survey Variables
		sumtable(dt_long,
			vars=c('Z_pop_assign','age','gender','educ_high','catholic','rank_spo','years_served','years_served_mun','lives_in_sorsogon','lives_where_works','top4match_pnplgu','officer_attitude_idx','officer_attitude_idx_alt','trust_idx','trust_idx_alt','trust_safe','trust_intent','trust_info','empathy_idx','empathy_idx_alt','empathy_complain','empathy_report','accountability_idx','accountability_idx_alt','account_care','h2_punish','h2_report_self','h2_report_other','h3_punish','h3_report_self','h3_report_other','corrupt_idx','corrupt_idx_alt','h2_corrupt_self','h2_corrupt_other','h3_corrupt_self','h3_corrupt_other'),
			labels=c('Treated','Age','Male','Completed College','Catholic','Rank (Senior)','Years Served','Years Served in Mun','Lives in Province','Lives in Municipality','Share Concerns','Officer Attitude Index','Officer Attitude Idx (alt)','Trust Index','Trust Index (alt)','Trust (Safe)','Trust (Intent)','Trust (Info)','Empathy Index','Empathy Index (alt)','Empathy (Complaints)','Empathy (Reports)','Accountability Index','Accountability Idx (alt)','Accountability (Care)','Accountability (Hyp1 - Punish)','Accountability (Hyp1 - Report Self)','Accountability (Hyp1 - Report Other)','Accountability (Hyp2 - Punish)','Accountability (Hyp2 - Report Self)','Accountability (Hyp2 - Report Other)','Corruption Idx','Corruption Idx (alt)','Corruption (Hyp1 - Self)','Corruption (Hyp1 - Other','Corruption (Hyp2 - Self)','Corruption (Hyp2 - Other'),			
			out='latex',file='3_Tables/0_raw/TabA1_SummaryStatsSvy_raw.tex',
			summ=c('notNA(x)','mean(x)','sd(x)'),group='midline',
			summ.names=c('N','Mean','SD'),digits=2,title='Summary Statistics - Survey'
			)

		# Municipal Variables
		mun.attr<-dt_long[!duplicated(station),c('station','c10_population','c10_hhsize','c10_education','c10_catholic','c10_home_good','npa_2015_presence')]
		sumtable(mun.attr,
			vars=c('c10_population','c10_hhsize','c10_education','c10_catholic','c10_home_good','npa_2015_presence'),
			labels=c('Population (mun)','Household Size (mun)','Education (mun)','Percent Catholic (mun)','Home Quality (mun)','Avg. NPA Presence 2015 (mun)'),			
			out='latex',file='3_Tables/0_raw/TabA2_SummaryStatsMun_raw.tex',
			summ=c('notNA(x)','mean(x)','sd(x)','min(x)','max(x)'),
			summ.names=c('N','Mean','SD','Min','Max'),digits=2,title='Summary Statistics - Municipalities'
			)

	# =======================================================
	# Descriptive - Balance Table
	# =======================================================	
		# Prep		
		balance <-
		  c("age", "gender", "educ_high",'catholic','lives_in_sorsogon','commute_m','years_served','top4match_pnplgu_m','officer_attitude_idx_m','trust_idx_m','empathy_idx_m','accountability_idx_m','corrupt_idx_m') %>%
		  map(~ formula(paste0('Z_pop_assign ~ ',.,'+ rank_m_spo'))) %>%
		  map(~ lm_robust(., data=dt[!is.na(sid_m)],clusters=station_m, fixed_effects= ~ station_m, se_type="stata")) %>%
		  map_df(tidy)
		balance<-as.data.table(balance)
		balance<-balance[term!='rank_m_spo',c('term','estimate','p.value')]
		balance$term<-c("Age", "Male", "Education",'Catholic','Lives in Sorsogon','Commute','Years Served','Share Citizen Concerns','Officer Attitude Index','Trust Index','Empathy Index','Accountability Index','Corruption Index')

		# Export Table
		kable(balance, digits =3, col.names=c('Variable','Estimate','p-value'),align = "lrr",format='latex',position='h') %>%
		save_kable('3_Tables/0_raw/TabA3_BalanceTable_raw.tex')

	# =======================================================
	# Descriptive - Attrition
	# =======================================================	
		# --- Tests for differential attrition (covariates) --- # 
	    	# Order to generate correct comparison groups		
	    	dt[,attrite_cat:=factor(attrite_cat,levels=c('Baseline Only','Endline Only','Both Surveys'))]	

			# Run Regressions w/ the covariates as DVs    			
			attrition1 <-
			  c('gender','age','educ_high','catholic','years_served','lives_in_sorsogon','npa_2015_presence','civ_crime_victim_idx_admin') %>%
			  map(~ formula(paste0(.,'~ attrite_cat '))) %>% map(~ lm_robust(., data=dt,clusters=station_m, se_type="stata")) %>% map_df(tidy)			

		# --- Tests for differential attrition (outcomes) --- # 
	    	# Order to generate correct comparison groups		
	    	dt[,attrite_cat:=factor(attrite_cat,levels=c('Both Surveys','Baseline Only','Endline Only'))]	

			# Run Regressions w/ outcome indices as DVs	
			attrition2 <-
			  c('officer_attitude_idx_m_attrite','trust_idx_m_attrite','empathy_idx_m_attrite','accountability_idx_m_attrite','corrupt_idx_m_attrite','officer_attitude_idx_e','trust_idx_e','empathy_idx_e','accountability_idx_e','corrupt_idx_e') %>%
			  map(~ formula(paste0(.,'~ attrite_cat'))) %>% map(~ lm_robust(., data=dt[Z_pop_assign==0], clusters=station_m,se_type="stata")) %>% map_df(tidy)			  

		# --- Export Table --- #
			# Prep
			attrition<-rbind(as.data.table(attrition1),as.data.table(attrition2))			  
			attrition<-attrition[term %!in% c('rank_m_spo','(Intercept)','Z_pop_assign'),c('outcome','term','estimate','p.value')]
			attrition[,p.value:=round(p.value,2)]; attrition[,estimate:=round(estimate,3)]
			attrition$outcome<-c('Male', "Male","Age",'Age', "College Educated",'College Educated','Catholic','Years Served','Years Served','Lives in Sorsogon','Lives in Sorsogon','NPA Presence 2015','NPA Presence 2015','Crime Rate','Crime Rate','Attitude Idx (baseline)','Trust Idx (baseline)','Empathy Idx (baseline)','Account Idx (baseline)','Corrupt Idx (baseline)','Attitude Idx (endline)','Trust Idx (endline)','Empathy Idx (endline)','Account Idx (endline)','Corrupt Idx (endline)')
			attrition[,term:=gsub('attrite_cat','',term)]

			# Export
			kable(attrition, digits =3, col.names=c('Variable','Attrition Category','Estimate','p-value'),align = "llrr",format='latex',position='h') %>%
			save_kable('3_Tables/0_raw/TabA4_Attrition_raw.tex')			

	# =======================================================
	# Descriptive - Balance of assignments on embeddedness
	# =======================================================			
		# Run regressions w/ the municipality-level characteristics as the DVs
		b1<-lm_robust(c10_education ~ lives_in_sorsogon, data=dt,clusters=station_m,se_type='stata') %>% tidy 
		b2<-lm_robust(c10_population ~ lives_in_sorsogon, data=dt,clusters=station_m,se_type='stata') %>% tidy 
		b3<-lm_robust(c10_home_good ~ lives_in_sorsogon, data=dt,clusters=station_m,se_type='stata') %>% tidy 
		b4<-lm_robust(c10_hhsize ~ lives_in_sorsogon, data=dt,clusters=station_m,se_type='stata') %>% tidy 
		b5<-lm_robust(c10_catholic ~ lives_in_sorsogon, data=dt,clusters=station_m,se_type='stata') %>% tidy 
		b6<-lm_robust(civ_crime_victim_idx_admin ~ lives_in_sorsogon, data=dt,clusters=station_m,se_type='stata') %>% tidy 
		b7<-lm_robust(npa_2015_presence ~ lives_in_sorsogon, data=dt,clusters=station_m,se_type='stata') %>% tidy 
		balance_embed<-as.data.table(rbind(b1,b2,b3,b4,b5,b6,b7))

		# Prep
		balance_embed<-balance_embed[term=='lives_in_sorsogon',c('outcome','estimate','p.value')]
		balance_embed[,p.value:=round(p.value,2)]; balance_embed[,estimate:=round(estimate,3)]
		balance_embed$outcome<-c('Municipality Education','Municipality Population','Municipality Home Quality','Municipality Household Size','Municipality Percent Catholic','Municipality Crime Rate','Municipality NPA Presence')

		# Export Table
		kable(balance_embed, digits =3, col.names=c('Outcome Variable','Estimate','p-value'),align = "lrr",format='latex',position='h') %>%
		save_kable('3_Tables/0_raw/TabA11_BalanceEmbed_raw.tex')		

	# =======================================================
	# Models - Main Average Treatment Effects
	# =======================================================	
		# Table
		texreg(
		  l=list(m_top4_panel,m_top4_cross,m_attitude_panel,m_attitude_cross,m_trust_panel,m_trust_cross,m_empathy_panel,m_empathy_cross,m_account_panel,m_account_cross,m_corrupt_panel,m_corrupt_cross),
		  file = '3_Tables/0_raw/TabA5_MainResults_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Share Concerns" = 1:2,' Attitude Idx'=3:4,"Trust" = 5:6,'Empathy'=7:8,'Accountability'=9:10,'Corruption'=11:12),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)','(8)','(9)','(10)','(11)','(12)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','top4match_pnplgu_m'='Base Outcome','trust_idx_m'='Base Outcome','empathy_idx_m'='Base Outcome','accountability_idx_m'='Base Outcome','officer_attitude_idx_m'='Base Outcome','corrupt_idx_m'='Base Outcome'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Main Treatment Effects", caption.above = TRUE, label = "tab:MainResults", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",
		  custom.gof.rows = list('Model'=c('Panel','Cross','Panel','Cross','Panel','Cross','Panel','Cross','Panel','Cross','Panel','Cross')),
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Panel models include a dummy variable for missing baseline. \n\\item All models include station (municipality) and rank fixed effects for the blocking variables. \n\\item Standard errors are clustered at the station-level.')

	# =======================================================
	# Models - Heterogeneous Effects: Embeddedness
	# =======================================================	
		# Without Controls
		texreg(
		  l=list(m_top4_embed,m_attitude_embed,m_trust_embed,m_empathy_embed,m_accountability_embed,m_corrupt_embed,m_attitude_alt_embed),
		  file = '3_Tables/0_raw/TabA6_HetEffects_Embed_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Share" = 1,"Attitude" = 2,"X" = 3,"X" = 4,"X" =5,"X"=6,"Attitude"=7),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','lives_in_sorsogon'='Lives in Province','Z_pop_assign:lives_in_sorsogon'='Treat*InProvince','top4match_pnplgu_m'='Base Outcome','trust_idx_m'='Base Outcome','empathy_idx_m'='Base Outcome','accountability_idx_m'='Base Outcome','officer_attitude_idx_m'='Base Outcome','corrupt_idx_m'='Base Outcome','officer_attitude_idx_alt_m'='Base Outcome'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Heterogeneous Treatment Effects - Embeddedness", caption.above = TRUE, label = "tab:HetEffects_Embed", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Models in this table use the panel specification denoted in equation 1 and include a dummy variable for missing baseline. All models also include station (municipality) and rank fixed effects for the blocking variables. Standard errors are clustered at the station-level.'	
		  )

		# With Controls
		texreg(
		  l=list(m_attitude_embed_ctrl1,m_attitude_embed_ctrl2,m_attitude_embed_ctrl3,m_trust_embed_ctrl3,m_empathy_embed_ctrl3,m_account_embed_ctrl3,m_corrupt_embed_ctrl3),
		  file = '3_Tables/0_raw/TabA9_HetEffects_EmbedwControl_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Attitude Index" = 1:3,"Trust"=4,"Empathy"=5,"Account"=6,"Corrupt"=7),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','lives_in_sorsogon'='Lives in Province','Z_pop_assign:lives_in_sorsogon'='Treated*InProvince','top4match_pnplgu_m'='Base Outcome','trust_idx_m'='Base Outcome','empathy_idx_m'='Base Outcome','accountability_idx_m'='Base Outcome','officer_attitude_idx_m'='Base Outcome','corrupt_idx_m'='Base Outcome','rank_m_spo'='Rank (officer)','educ_high'='Education (officer)','catholic'='Catholic (officer)','age'='Age (officer)','gender'='Gender (officer)','years_served_mun'='Years Served (officer)','commute_m'='Commute (officer)','psm_idx_combined_m'='PSM (officer)','c10_education'='Education (mun)','c10_population'='Population (mun)','c10_homegood'='Home Quality (mun)','c10_hhsize'='Household Size (mun)','c10_catholic'='Percent Catholic (mun)','civ_crime_victim_idx_admin'='Crime Rate (mun)','npa_2015_presence'='NPA Presence 2015 (mun)','Z_pop_assign:rank_m_spo'='Treat*Rank','Z_pop_assign:educ_high'='Treat*Education','Z_pop_assign:catholic'='Treat*Catholic','Z_pop_assign:age'='Treat*Age','Z_pop_assign:gender'='Treat*Gender','Z_pop_assign:years_served_mun'='Treat*YearsServed','Z_pop_assign:commute_m'='Treat*Commute','Z_pop_assign:psm_idx_combined_m'='Treat*PSM'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Heterogeneous Treatment Effects - Embeddedness (with Controls)", caption.above = TRUE, label = "tab:HetEffects_Embed_wControl", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",center=TRUE,
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Models in this table use the panel specification denoted in equation 1 and include a dummy variable for missing baseline. Standard errors are clustered at the station-level.'
		  )

		# Alternate Embeddedness Measures
		texreg(
		  l=list(m_top4_embedmun,m_attitude_embedmun,m_trust_embedmun,m_empathy_embedmun,m_accountability_embedmun,m_corrupt_embedmun,m_attitude_alt_embedmun),
		  file = '3_Tables/0_raw/TabA14_HetEffects_EmbedMun_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
			  custom.header = list("Share" = 1,"Attitude" = 2,"X" = 3,"X" = 4,"X" =5,"X"=6,"Attitude"=7),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','lives_where_works'='Lives in Municipality','Z_pop_assign:lives_where_works'='Treated*InMunicipality','top4match_pnplgu_m'='Base Outcome','trust_idx_m'='Base Outcome','empathy_idx_m'='Base Outcome','accountability_idx_m'='Base Outcome','officer_attitude_idx_m'='Base Outcome','officer_attitude_idx_alt_m'='Base Outcome','corrupt_idx_m'='Base Outcome','rank_m_spo'='Rank (officer)','educ_high'='Education (officer)','catholic'='Catholic (officer)','age'='Age (officer)','gender'='Gender (officer)','years_served_mun'='Years Served (officer)','commute_m'='Commute (officer)','psm_idx_combined_m'='PSM (officer)'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Heterogeneous Treatment Effects - Embeddedness (Lives In Municipality)", caption.above = TRUE, label = "tab:HetEffects_EmbedMun", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",center=TRUE,
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Models in this table use the panel specification denoted in equation 1 and include a dummy variable for missing baseline. Models include station and rank fixed effects. Standard errors are clustered at the station-level.'
		  )

	# =======================================================
	# Models - Heterogeneous Effects: NPA
	# =======================================================	
		# Main Outcomes 
		texreg(
		  l=list(m_top4_npa,m_attitude_npa,m_trust_npa,m_empathy_npa,m_account_npa,m_corrupt_npa,m_trust_alt_npa),
		  file = '3_Tables/0_raw/TabA8_HetEffects_NPA_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Share" = 1,"Attitude" = 2,"X" = 3,"X" = 4,'X'=5,'X'=6,'Trust'=7),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','npa_2015_presence'='NPA Presence 2015','Z_pop_assign:npa_2015_presence'='Treated*NPA','top4match_pnplgu_m'='Base Outcome','trust_idx_m'='Base Outcome','trust_idx_m'='Base Outcome','empathy_idx_m'='Base Outcome','accountability_idx_m'='Base Outcome','officer_attitude_idx_m'='Base Outcome','corrupt_idx_m'='Base Outcome'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Heterogeneous Treatment Effects - Rebel Presence", caption.above = TRUE, label = "tab:HetEffects_NPA", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",center=TRUE,
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Models in this table use the panel specification denoted in equation 1 and include a dummy variable for missing baseline. All models also include rank fixed effects for the blocking variables, though they do not include station fixed effects because NPA Presence is measured at the municipality (station) level. Standard errors are clustered at the station-level.'
		  )

		# Trust Sub-Indices and NPA Sympathy
		texreg(
		  l=list(m_trust_npa,m_tsafe_npa,m_tintent_npa,m_tinfo_npa,m_sympathy_npa),
		  file = '3_Tables/0_raw/TabA7_HetEffects_NPAsub_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Trust" = 1,"Trust" = 2,"Trust" = 3,"Trust" = 4,'Perceived'=5),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','npa_2015_presence'='NPA Presence 2015','Z_pop_assign:npa_2015_presence'='Treated*NPA','trust_idx_m'='Base Outcome','trust_idx_m'='Base Outcome','trust_safe_m'='Base Outcome','trust_intent_m'='Base Outcome','trust_info_m'='Base Outcome','npa_sympathy_m'='Base Outcome'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Heterogeneous Treatment Effects - Rebel Presence (Trust Outcomes)", caption.above = TRUE, label = "tab:HetEffects_NPA", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",center=TRUE,
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Models in this table use the panel specification denoted in equation 1 and include a dummy variable for missing baseline. All models also include rank fixed effects for the blocking variables, though they do not include station fixed effects because NPA Presence is measured at the municipality (station) level. Standard errors are clustered at the station-level.'
		  )

		# Trust Sub-Indices and NPA Sympathy (With Controls)
		texreg(
		  l=list(m_trust_npa_ctrl1,m_trust_npa_ctrl2,m_trust_npa_ctrl3,m_tsafe_npa_ctrl2,m_tintent_npa_ctrl2,m_tinfo_npa_ctrl2,m_trust_alt_npa_ctrl2,m_sympathy_npa_ctrl2,m_sympathy_npa_ctrl1),
		  file = '3_Tables/0_raw/TabA10_HetEffects_NPAwControl_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Trust" = 1:3,"Trust" = 4,"Trust" = 5,"Trust" = 6,"Trust" = 7,'Perceived'=8:9),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)','(8)','(9)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','npa_2015_presence'='NPA Presence 2015','Z_pop_assign:npa_2015_presence'='Treated*NPA','trust_idx_m'='Base Outcome','trust_idx_alt_m'='Base Outcome','trust_safe_m'='Base Outcome','trust_intent_m'='Base Outcome','trust_info_m'='Base Outcome','npa_sympathy_m'='Base Outcome','c10_education'='Education (mun)','c10_population'='Population (mun)','c10_homegood'='Home Quality (mun)','c10_hhsize'='Household Size (mun)','c10_catholic'='Percent Catholic (mun)','civ_crime_victim_idx_admin'='Crime Rate (mun)','rank_m_spo'='Rank (officer)','educ_high'='Education (officer)','catholic'='Catholic (officer)','age'='Age (officer)','gender'='Gender (officer)','years_served_mun'='Years Served (officer)','commute_m'='Commute (officer)','psm_idx_combined_m'='PSM (officer)'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Heterogeneous Treatment Effects - Rebel Presence (with Controls)", caption.above = TRUE, label = "tab:HetEffects_NPA_wControl", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",center=TRUE,
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Models in this table use the panel specification denoted in equation 1 and include a dummy variable for missing baseline. Models do not include station fixed effects because NPA Presence is measured at the municipality (station) level. Standard errors are clustered at the station-level.'
		  )

	# =======================================================
	# Models - Robustness: Alternative Measures
	# =======================================================	
		# Alternative Attitude Indices -- Main Models
		texreg(
		  l=list(m_attitude_alt_panel,m_attitude_alt_cross,m_trust_alt_panel,m_trust_alt_cross,m_empathy_alt_panel,m_empathy_alt_cross,m_account_alt_panel,m_account_alt_cross,m_corrupt_alt_panel,m_corrupt_alt_cross),
		  file = '3_Tables/0_raw/TabA12_MainResultsAlt_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Attitude Idx" = 1:2,'Trust'=3:4,'Empathy'=5:6,'Accountability'=7:8,'Corruption'=9:10),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)','(8)','(9)','(10)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','trust_idx_alt_m'='Base Outcome','empathy_idx_alt_m'='Base Outcome','accountability_idx_alt_m'='Base Outcome','officer_attitude_idx_alt_m'='Base Outcome','corrupt_idx_alt_m'='Base Outcome'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Treatment Effects - Alternate Indices", caption.above = TRUE, label = "tab:MainResults_alt", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",
		  custom.gof.rows = list('Model'=c('Panel','Cross','Panel','Cross','Panel','Cross','Panel','Cross','Panel','Cross')),
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Panel models include a dummy variable for missing baseline. \n\\item All models include station (municipality) and rank fixed effects for the blocking variables. \n\\item Standard errors are clustered at the station-level.')

		# Alternative Knowledge Measures
		texreg(
		  l=list(m_top4_panel,m_top4_cross,m_top3_panel,m_top3_cross,m_top4know_panel,m_top4know_cross,m_top4civ_panel,m_top4civ_cross),
		  file = '3_Tables/0_raw/TabA13_KnowledgeAlt_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Top 4 LGU/own" = 1:2,"Top 3 LGU/own" = 3:4,'Top 4 LGU/their'=5:6,'Top 4 CIV/own'=7:8),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)','(8)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','top4match_pnplgu_m'='Base Outcome','top3match_pnplgu_m'='Base Outcome','top4match_civlgu_m'='Base Outcome','top3match_pnpciv_m'='Base Outcome'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Treatment Effects: Alternative Knowledge Measures", caption.above = TRUE, label = "tab:KnowledgeAlt", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",
		  custom.gof.rows = list('Model'=c('Panel','Cross','Panel','Cross','Panel','Cross','Panel','Cross')),
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Panel models include a dummy variable for missing baseline. \n\\item All models include station (municipality) and rank fixed effects for the blocking variables \n\\item Models labeled with "Top 4" indicate that the outcome was measured as the match of officer responses with the top 4 issues mentioned by LGU/civilian respondents in the same municipality, while models labeled "Top 3" indicate that we matched with only the top 3 issues mentioned. Models labeled "LGU" indicate that we matched officer responses with the aggregate responses of the LGU community leaders, while models labeled "CIV" indicate that we matched with the aggregate responses of ordinary civilians. Models labeled "own" indicate that we matched using the issues that officers themselves thought were important, while models labeled "their" indicate that we matched using the issues that officers thought that LGU/civilian respondents thought were important.')

	# =======================================================
	# Models - Spillovers
	# =======================================================
		# Station-Level Moderators
		texreg(
		  l=list(m_top4_statsize,m_attitude_statsize,m_trust_statsize,m_empathy_statsize,m_accountability_statsize,m_corruption_statsize,m_top4_statdensity,m_attitude_statdensity,m_trust_statdensity,m_empathy_statdensity,m_accountability_statdensity,m_corruption_statdensity),
		  file = '3_Tables/0_raw/TabA24_Spillovers_Station_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list('Share'=1,'Attitude'=2,'X'=3,'X'=4,'X'=5,'X'=6,'Share'=7,'Attitude'=8,'X'=9,'X'=10,'X'=11,'X'=12),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)','(8)','(9)','(10)','(11)','(12)'),
		  custom.coef.map = list('Z_pop_assign'='Treated','station_size'='StationSize','Z_pop_assign:station_size'='Treated*Size','density_station'='StationDensity','Z_pop_assign:density_station'='Treated*Density','top4match_pnplgu_m'='Base Outcome','trust_idx_alt_m'='Base Outcome','trust_idx_m'='Base Outcome','empathy_idx_m'='Base Outcome','accountability_idx_m'='Base Outcome','officer_attitude_idx_m'='Base Outcome','corrupt_idx_m'='Base Outcome','c10_education'='Education','c10_population'='Population','c10_homegood'='Home Quality','c10_hhsize'='HH Size','c10_catholic'='Catholic','civ_crime_victim_idx_admin'='Crime Rate'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 3,caption = "Spillover Analysis - Station-Level", caption.above = TRUE, label = "tab:spillovers_station", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item Models in this table use the panel specification denoted in equation 1 and include a dummy variable for missing baseline. \n\\item All models include rank fixed effects, but do not include station fixed effects because the moderating variables are measured at the station level. \n\\item Standard errors are clustered at the station-level.'
		  )

		# Individual-Level Moderators
		texreg(
		  l=list(m_knowledge_spillany,m_knowledge_spillfam,m_knowledge_spillbeat,m_attitude_spillany,m_attitude_spillfam,m_attitude_spillbeat,m_trust_spillany,m_empathy_spillany,m_accountability_spillany,m_corrupt_spillany),
		  file = '3_Tables/0_raw/TabA23_Spillovers_Individual_raw.tex',
		  stars = c(0.01, 0.05, 0.10),
		  custom.header = list("Shared Concerns" = 1:3,'Attitude Index'=4:6,'Trust'=7,'Empathy'=8,'Account'=9,'Corrupt'=10),
		  custom.model.names = c('(1)','(2)','(3)','(4)','(5)','(6)','(7)','(8)','(9)','(10)'),
		  custom.coef.map = list('spill_any'='Exposure to treated (Any)','spill_family_direct'='Exposure to treated (Family)','spill_beat_any'='Exposure to treated (Beat)'),
		  include.adjrs = FALSE,include.rmse=FALSE, include.ci = FALSE,digits = 2,caption = "Spillover Analysis - Individual-Level", caption.above = TRUE, label = "tab:spillovers_ind", booktabs = TRUE, sideways = FALSE, no.margin = TRUE, use.packages=FALSE, threeparttable=TRUE, dcolumn=TRUE, custom.gof.names = c(NA,'Observations','Clusters'),float.pos = "h",
		  custom.note = '\\item \\sym{*} \\(p<.10\\), \\sym{**} \\(p<.05\\), \\sym{***} \\(p<.01\\)} \n\\item All models include station (municipality) and rank fixed effects for the blocking variables \n\\item Regressions are subset to officers in the control group. The independent variables indicate whether the officer had a family member in the treatment group (Family), a beat patrol partner in the treatment group (Beat), or either a family member or beat patrol partner in the treatment group (Any).'
		  )




















